-Aleksander Nistad Sekkelsten
import Pkg; Pkg.activate(".julia/environments/v1.10/Project.toml")
Pkg.add(["SCS", "Convex","PlotThemes","Images"])
using Flux, Plots, Random, NBInclude, LinearAlgebra, Statistics, Distributions, Convex, SCS, Images
using Plots.PlotMeasures
using StatsBase: sample
theme(:dao)
Activating project at `~/.julia/environments/v1.10` Updating registry at `~/.julia/registries/General.toml` Resolving package versions... No Changes to `~/.julia/environments/v1.10/Project.toml` No Changes to `~/.julia/environments/v1.10/Manifest.toml`
@nbinclude("/Users/aleksandersekkelsten/Downloads/Project 1 function.ipynb")
CrossValidation
?CrossValidation
search: CrossValidation
CrossValidation(X, Y, Z, k_folds::Int64, degrees::Int64, Method; lambda=0)
Performs k-fold cross-validation to evaluate a regression model on the provided dataset. The function splits the data into k_folds, fits a model on the training data, and evaluates it on the test data for each fold, returning the mean error across all folds.
X::Array{Float64}: The x-coordinates of the input data (for grid).
Y::Array{Float64}: The y-coordinates of the input data (for grid).
Z::Array{Float64}: The z-values (target) of the dataset.
k_folds::Int64: The number of folds to use in cross-validation. The data is split into k_folds parts, where each part is used once as the test set.
degrees::Int64: The degree of the polynomial used to construct the design matrix.
Method::Symbol: The regression method to use. Can be:
:OLS for Ordinary Least Squares.:Ridge for Ridge regression.:Lasso for Lasso regression.lambda::Float64 (optional): The regularization parameter for Ridge or Lasso regression. Default is 0 (no regularization).
mean_error::Float64: The mean error (MSE) across all the k-folds of cross-validation.Data Processing: The input data X, Y, and Z are flattened and reshaped to prepare for the regression model. A design matrix is constructed using the DataProjector function, mapping the inputs into a higher-dimensional polynomial space defined by degrees.
Shuffling and Splitting: The data is shuffled randomly to remove any bias from its original ordering. The shuffled data is split into k_folds parts for cross-validation.
Cross-Validation Loop: For each fold:
Error Calculation: The mean squared error (MSE) is computed for each fold, and the average error across all folds is returned as the overall performance metric.
meanerror = CrossValidation(X, Y, Z, k_folds, degrees, Method, lambda=0.1)
#X and Y range
X = range(0,1,100) |> collect
Y = range(0,1,100) |> collect
Random.seed!(123)
ϵ=Normal(0, 0.02) # Noise
Z = Franke.(X',Y) .+ rand(ϵ,(100,100)) #Finding Z values
#Display Frankes funciton
surface(X,Y,Z, size=(1000,700),title="Franke's function")
DM,Z_flat = DataProjector(X,Y, Z, 5,train=false)
xTrain, yTrain, zTrain, xTest, yTest, zTest = SplitData2D(X,Y,Z,0.8)
D_M_train = DataProjector(xTrain, yTrain, zTrain, 5)
D_M_test = DataProjector(xTest, yTest, zTest, 5)
betaVals = Regression(D_M_train,zTrain,"OLS")
est = DM*betaVals
r2 = round(R2(D_M_test, zTest,betaVals),digits=4)
mse = round(MSE(D_M_test,zTest,betaVals),digits=4)
surface(X,Y,Z, layout=(1,2), cbar=false, c=:grays,size=(900,400),title="Ground Truth")
surface!(subplot=2,X,Y,reshape(est, (100,100)), cbar=false, title="Regression, MSE: $(mse) R2 = $(r2)")
DM,Z_flat = DataProjector(X,Y, Z, 5,train=false,normalize=false)
([1.0 0.0 … 0.0 0.0; 1.0 0.010101010101010102 … 0.0 0.0; … ; 1.0 0.98989898989899 … 0.98989898989899 1.0; 1.0 1.0 … 1.0 1.0], [0.7825863498542225, 0.7408899906204466, 0.7377043500767736, 0.7474868504560255, 0.7571751612607839, 0.7511803004298346, 0.732863392429973, 0.7084405670141984, 0.730712831513705, 0.7202400702570643 … -0.06604138053076845, -0.012511207278854264, -0.02760643785794465, -0.0058408863780550465, -0.0029773028561855567, -0.012356191438078104, -0.016299631557830572, 0.0025389500852499142, -0.01584425417806084, 0.02121899090848535])
lambdaVals = [1e-4,1e-1,1]
#Solve β for OLS and find estimate
betaValsOLS = Regression(D_M_train,zTrain,"OLS")
estOLS = DM*betaValsOLS
#R2 and MSE value for OLS
r2OLS = round(R2(D_M_test, zTest,betaValsOLS),digits=4)
mseOLS = round(MSE(D_M_test,zTest,betaValsOLS),digits=4)
println("OLS | MSE = ", mseOLS, " R2 value = ", r2OLS, "\n")
surf = surface(X,Y,Z, layout=(4,2), cbar=false, c=:grays,size=(500,1000),title="Ground Truth")
surface!(surf,subplot=2,X,Y,reshape(estOLS, (100,100)), cbar=false, title="OLS Regression ")
for (i,lambda) in enumerate(lambdaVals) #Looping over various λ values: lambdaVals = [1e-4,1e-2,1e-1]
#Solve β for Lasso and Ridge
betaLasso = Regression(D_M_train, zTrain, "Lasso", λ=lambda)
betaRidge = Regression(D_M_train, zTrain, "Ridge", λ=lambda)
#Calculate MSE
MSELasso = round(MSE(D_M_test, zTest, betaLasso, λ=0, γ=lambda),digits=4)
MSERidge = round(MSE(D_M_test, zTest, betaRidge, λ=lambda, γ=0),digits=4)
#Calculate R2 value
r2Lasso = round(R2(D_M_test, zTest, betaLasso),digits=4)
r2Ridge = round(R2(D_M_test, zTest, betaRidge),digits=4)
println("Lasso | λ = ", lambda," MSE: ", MSELasso, " R2 value = ", r2Lasso)
println("Ridge | λ = ", lambda," MSE: ", MSERidge, " R2 value = ", r2Ridge, "\n")
#Find estimated answer
estLasso = DM*betaLasso
estRidge = DM*betaRidge
#Display result
surface!(surf,subplot=2*(i+1),X,Y,reshape(estLasso, (100,100)), cbar=false, title="Lasso Regression | λ= $(lambda)")
surface!(surf,subplot=(2*(i+1)-1),X,Y,reshape(estRidge, (100,100)), cbar=false, title="Ridge Regression | λ= $(lambda)")
end
OLS | MSE = 0.0025 R2 value = 0.9704 Lasso | λ = 0.0001 MSE: 0.0414 R2 value = 0.9704 Ridge | λ = 0.0001 MSE: 1.0328 R2 value = 0.9704 Lasso | λ = 0.1 MSE: 10.7632 R2 value = 0.9559 Ridge | λ = 0.1 MSE: 11.133 R2 value = 0.9291 Lasso | λ = 1.0 MSE: 22.8137 R2 value = 0.9122 Ridge | λ = 1.0 MSE: 12.4153 R2 value = 0.8902
display(surf)
This section is rather large, we first compare Test/Train Errors, then introduce Bootstrap resampling and analyze Error, Bias and Variance. Then we look at Cross-Validation. (Task e-f)
degrees = range(1,13) |> collect #Range of polynomial degrees
Test_Errors = zeros(length(degrees))
Train_Errors = zeros(length(degrees))
for (i,degree) in enumerate(degrees)
# ---- Initializing Data ----
xTrain, yTrain, zTrain, xTest, yTest, zTest = SplitData2D(X,Y,Z,0.7)
D_M_train = DataProjector(xTrain, yTrain, zTrain, degree)
D_M_test = DataProjector(xTest, yTest, zTest, degree)
# ---- Solving OLS ----
betaValsOLS = Regression(D_M_train,zTrain,"OLS")
# ---- Updating Errors ----
Test_Errors[i] = MSE(D_M_test,zTest,betaValsOLS)
Train_Errors[i] = MSE(D_M_train,zTrain,betaValsOLS)
end
plot(degrees, log.(Test_Errors), label="Test Errors", xlabel="Model Complexity", ylabel="Error", title = "OLS Regression")
plot!(degrees, log.(Train_Errors), label="Train Errors")
lambda = 1e-3 #λ for Ridge Regression
degrees = range(1,25) |> collect #Range of polynomial degrees
Test_Errors = zeros(length(degrees))
Train_Errors = zeros(length(degrees))
Bias = zeros(length(degrees))
Variance = zeros(length(degrees))
for (i,degree) in enumerate(degrees)
# ---- Initializing Data ----
xTrain, yTrain, zTrain, xTest, yTest, zTest = SplitData2D(X,Y,Z,0.7)
D_M_train = DataProjector(xTrain, yTrain, zTrain, degree)
D_M_test = DataProjector(xTest, yTest, zTest, degree)
# ---- Solving Ridge ----
betaVals = Regression(D_M_train,zTrain,"Ridge", λ=lambda)
# ---- Updating Errors ----
Test_Errors[i] = MSE(D_M_test,zTest,betaVals, λ = lambda,γ=0)
Train_Errors[i] = MSE(D_M_train,zTrain,betaVals, λ = lambda, γ=0)
end
plot(degrees, Test_Errors, label="Test Errors", xlabel="Model Complexity", ylabel="Error", title="Ridge Regression")
plot!(degrees, Train_Errors, label="Train Errors")
lambda = 1e-3 #λ for Ridge Regression
degreeList = range(1,25) |> collect #Range of polynomial degrees
BiasOLS = zeros(length(degrees))
VarianceOLS = zeros(length(degrees))
BiasRidge = zeros(length(degrees))
VarianceRidge = zeros(length(degrees))
BiasLasso = zeros(length(degrees))
VarianceLasso = zeros(length(degrees))
for (i, degree) in enumerate(degreeList)
# ---- Initializing Data ----
xTrain, yTrain, zTrain, xTest, yTest, zTest = SplitData2D(X,Y,Z,0.7)
D_M_train = DataProjector(xTrain, yTrain, zTrain, degree)
D_M_test = DataProjector(xTest, yTest, zTest, degree)
# ---- Solving OLS, Ridge, Lasso ----
betaValsOLS = Regression(D_M_train,zTrain,"OLS")
betaValsRidge = Regression(D_M_train,zTrain,"Ridge", λ=lambda)
betaValsLasso = Regression(D_M_train,zTrain,"Lasso", λ=1e-2)
if length(betaValsLasso)==1
betaValsLasso = Regression(D_M_train,zTrain,"Lasso", λ=1e-3)
end
# ---- Calculating Bias & Variance ----
BiasOLS[i] = mean((mean(D_M_test*betaValsOLS) .-zTest).^2)
VarianceOLS[i] = var(D_M_test*betaValsOLS)
BiasRidge[i] = mean((mean(D_M_test*betaValsRidge) .-zTest).^2)
VarianceRidge[i] = var(D_M_test*betaValsRidge)
BiasLasso[i] = mean((mean(D_M_test*betaValsLasso) .-zTest).^2)
VarianceLasso[i] = var(D_M_test*betaValsLasso)
end
plot(layout=(1,3), degrees, (BiasOLS), label="Bias", xlabel="Model Complexity", ylabel="Error", title = "OLS Regression",size=(1200,300))
plot!(degrees, (VarianceOLS), label="Variance")
plot!(subplot=2, degrees, (BiasRidge), label="Bias", xlabel="Model Complexity", ylabel="Error", title = "Ridge Regression")
plot!(subplot=2, degrees, (VarianceRidge), label="Variance")
plot!(subplot=3, degrees, (BiasLasso), label="Bias", xlabel="Model Complexity", ylabel="Error", title = "Lasso Regression")
plot!(subplot=3,degrees, (VarianceLasso), label="Variance")
degs = range(1,12) |> collect
numBootstraps=1000
errors = zeros(length(degs))
biases = zeros(length(degs))
variances = zeros(length(degs))
for (i,deg) in enumerate(degs)
# ---- Initializing Data ----
D_M_train = DataProjector(xTrain, yTrain, zTrain, deg)
D_M_test = DataProjector(xTest, yTest, zTest, deg)
# ---- Apply Bootstrap ----
bootstrapBetas = Bootstrap(xTrain, yTrain, zTrain, deg, numBootstraps,"OLS")
# ---- Analyzing Results ----
est_test = zeros(length(zTest),numBootstraps)
for i in 1:numBootstraps
est_test[:,i] = D_M_test*bootstrapBetas[:,i]
end
# ---- Update: Error, Bias and Variance ----
errors[i] = mean( mean((zTest .- est_test).^2, dims=1) )
biases[i] = mean( (zTest .- mean(est_test, dims=1)).^2 )
variances[i] = mean( var(est_test, dims=1) )
end
degs = range(1,12) |> collect
plot(degs, log.(errors), label="Errors w/ OLS")
plot!(degs, log.(biases), label="Biases^2 w/ OLS")
plot!(degs, log.(variances), label="Variances w/ OLS")
p = surface(X,Y,Z, layout=(2,2),title="Ground truth",cbar=false)
for i in 10:12 #Degrees 10,11 and 12
# ---- Initializing Data ----
DM,Z_flat = DataProjector(X,Y, Z, i,train=false)
xTrain, yTrain, zTrain, xTest, yTest, zTest = SplitData2D(X,Y,Z,0.7)
D_M_train = DataProjector(xTrain, yTrain, zTrain, i)
D_M_test = DataProjector(xTest, yTest, zTest, i)
# ---- Solve for OLS ----
betaVals = Regression(D_M_train,zTrain,"OLS")
estim = DM*betaVals
# ---- Display ----
surface!(p,subplot=(i-8),X,Y,reshape(estim, (100,100)),title="Degree = $(i)",cbar=false)
end
display(p)
degs = range(1,15) |> collect
Ridge_errors = zeros(length(degs))
Ridge_biases = zeros(length(degs))
Ridge_variances = zeros(length(degs))
for (i,deg) in enumerate(degs)
# ---- Initializing Data ----
D_M_train = DataProjector(xTrain, yTrain, zTrain, deg)
D_M_test = DataProjector(xTest, yTest, zTest, deg)
# ---- Apply Bootstrap ----
bootstrapBetas = Bootstrap(xTrain, yTrain, zTrain, deg, numBootstraps,"Ridge", λ=1e-3)
# ---- Analyzing Results ----
est_test = zeros(length(zTest),numBootstraps)
for i in 1:numBootstraps
est_test[:,i] = D_M_test*bootstrapBetas[:,i]
end
# ---- Update: Error, Bias and Variance ----
Ridge_errors[i] = mean( mean((zTest .- est_test).^2, dims=1) )
Ridge_biases[i] = mean( (zTest .- mean(est_test, dims=1)).^2 )
Ridge_variances[i] = mean( var(est_test, dims=1) )
end
Ridge = plot(degs, (Ridge_errors), label="Errors w/ Ridge")
plot!(Ridge,degs, (Ridge_biases.^2), label="Biases^2 w/ Ridge")
plot!(Ridge, degs, (Ridge_variances), label="Variances w/ Ridge")
degrees = [10, 12, 25]
p = surface(X,Y,Z, layout=(2,2),title="Ground truth",cbar=false)
for (i,degree) in enumerate(degrees) #Degrees 10,11 and 12
# ---- Initializing Data ----
DM,Z_flat = DataProjector(X,Y, Z, degree,train=false)
xTrain, yTrain, zTrain, xTest, yTest, zTest = SplitData2D(X,Y,Z,0.7)
D_M_train = DataProjector(xTrain, yTrain, zTrain, degree)
D_M_test = DataProjector(xTest, yTest, zTest, degree)
# ---- Solve for Ridge ----
betaVals = Regression(D_M_train,zTrain,"Ridge", λ=1e-3)
estim = DM*betaVals
# ---- Display ----
surface!(p,subplot=(i+1),X,Y,reshape(estim, (100,100)),title="Degree = $(degree)",cbar=false)
end
display(p)
degree=12
CVOLS = CrossValidation(X, Y, Z, 10, degree,"OLS")
xTrain, yTrain, zTrain, xTest, yTest, zTest = SplitData2D(X,Y,Z,0.7)
D_M_train = DataProjector(xTrain, yTrain, zTrain, degree)
D_M_test = DataProjector(xTest, yTest, zTest, degree)
bootstrapBetas = Bootstrap(xTrain, yTrain, zTrain, degree, 1000,"OLS")
BOLS = MSE(D_M_test,zTest,vec(mean(bootstrapBetas,dims=2)), λ=0, γ=0)
CVRidge = CrossValidation(X, Y, Z, 10, degree, "Ridge", lambda=1e-4)
CVLasso = CrossValidation(X, Y, Z, 10, degree, "Lasso", lambda=1e-4)
println(" ----- CrossValidation Scores ----- ", "\n")
println(" OLS : ", CVOLS, " B-Strap : ", BOLS)
println(" Ridge : ", round(CVRidge,digits=4))
println(" Lasso : ", round(CVLasso,digits=4))
----- CrossValidation Scores ----- OLS : 1.607353698836731 B-Strap : 0.09148002642683092 Ridge : 1.6612 Lasso : 0.1854
CVOLS = CrossValidation(X, Y, Z, 10, 5,"OLS")
xTrain, yTrain, zTrain, xTest, yTest, zTest = SplitData2D(X,Y,Z,0.7)
D_M_train = DataProjector(xTrain, yTrain, zTrain, 5)
D_M_test = DataProjector(xTest, yTest, zTest, 5)
bootstrapBetas = Bootstrap(xTrain, yTrain, zTrain, 5, 1000,"OLS")
BOLS = MSE(D_M_test,zTest,vec(mean(bootstrapBetas,dims=2)), λ=0, γ=0)
CVRidge = CrossValidation(X, Y, Z, 10, 5, "Ridge", lambda=1e-4)
CVLasso = CrossValidation(X, Y, Z, 10, 5, "Lasso", lambda=1e-4)
println(" ----- CrossValidation Scores ----- ", "\n")
println(" OLS : ", CVOLS, " B-Strap : ", BOLS)
println(" Ridge : ", round(CVRidge,digits=4))
println(" Lasso : ", round(CVLasso,digits=4))
----- CrossValidation Scores ----- OLS : 0.00251552069966107 B-Strap : 0.0024708075394363318 Ridge : 1.0342 Lasso : 0.0413
# ---- Load Image ----
img = load("/Users/aleksandersekkelsten/Downloads/SRTM_data_Norway_2.tif")
img_array = convert(Array{Float64}, img) #3601 x 1801
# ---- Defining X and Y range ----
y_img = range(1,3601) |> collect
x_img = range(1,1801) |> collect
1801-element Vector{Int64}:
1
2
3
4
5
6
7
8
9
10
11
12
13
⋮
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
heat = heatmap(img_array, color = :greys, cbar=false)
# ---- Selecting an area of the Image ----
(xmin,xmax) = (1250,1350)
(ymin,ymax) = (2350,2450)
img=img_array[ymin:ymax,xmin:xmax]
y_img = range(ymin,ymax) |> collect
x_img = range(xmin,xmax) |> collect
# ---- Displaying The Selected Area ----
heat = heatmap(layout = (1,2), img_array[1:end,1:end], color = :greys, cbar=false,size=(1000,400))
plot!(subplot=1,heat,[xmin,xmax], [ymin,ymin], ls=:dot,c=:red,alpha=0.8,label=false)
plot!(subplot=1,heat,[xmin,xmax], [ymax,ymax], ls=:dot,c=:red,alpha=0.8,label=false)
plot!(subplot=1,heat,[xmin,xmin], [ymin,ymax], ls=:dot,c=:red,alpha=0.8,label=false)
plot!(subplot=1,heat,[xmax,xmax], [ymin,ymax], ls=:dot,c=:red,alpha=0.8,label=false)
heatmap!(subplot=2,img_array[ymin:ymax,xmin:xmax], color = :greys, cbar=false)
display(heat)
display(heat)
# ---- Initializing Data ----
deg=110
DM_img,img_flat = DataProjector(x_img,y_img, img, deg,train=false,normalize=true)
xTrain, yTrain, zTrain, xTest, yTest, zTest = SplitData2D(x_img,y_img, img,0.8)
D_M_train, Train_stds, Train_mus = DataProjector(xTrain, yTrain, zTrain, deg,normalize=true)
D_M_test, Test_stds, Test_mus = DataProjector(xTest, yTest, zTest, deg,normalize=true)
# ---- Solving β for OLS, Ridge and Lasso ----
betasRidge = Regression(D_M_train,zTrain,"Ridge", λ=1e-2)
estRidge = DM_img*(betasRidge)
betasOLS = Regression(D_M_train,zTrain,"OLS")
estOLS = DM_img*(betasOLS)
betasLasso = Regression(D_M_train,zTrain,"Lasso", λ=1e-2)
estLasso = DM_img*(betasLasso)
10201-element Vector{Float64}:
0.006796960225891766
0.006363264221755655
0.006353511371169007
0.005808360733627452
0.006643020964298147
0.0069090654630146395
0.0068182488245664705
0.006132252228227034
0.005216454767041299
0.006039055328625918
0.006068875902987353
0.005673000443596979
0.005503106081437785
⋮
0.008311505886501378
0.009432032909462264
0.008819216331134548
0.00910466090439456
0.008819224240008952
0.009481123682981011
0.009481199318270335
0.00910733765578676
0.009045331892089612
0.00953627399808817
0.008829880293557683
0.00914050394036022
heatmap(layout=(1,3),reshape(estOLS, (ymax-ymin+1,xmax-xmin+1)), color = :greys,size=(2000,500),cbar=false,title="OLS",titlefontsize=25, top_margin=10mm)
heatmap!(subplot=2,reshape(estRidge, (ymax-ymin+1,xmax-xmin+1)), color = :greys,cbar=false,title="Ridge",titlefontsize=25, top_margin=10mm)
heatmap!(subplot=3,reshape(estLasso, (ymax-ymin+1,xmax-xmin+1)), color = :greys,cbar=false,title="Lasso",titlefontsize=25, top_margin=10mm)